Quantum Monte Carlo study of dilute neutron matter at finite temperatures 
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We report results of fully non-perturbative, Path Integral Monte Carlo (PIMC) calculations for 
dilute neutron matter. The neutron-neutron interaction in the s channel is parameterized by the 
scattering length and the effective range. We calculate the energy and the chemical potential as 
a function of temperature at the density p = 0.003 fm~ 3 . The critical temperature T c for the 
superfluid-normal phase transition is estimated from the finite size scaling of the condensate frac- 
tion. At low temperatures we extract the spectral weight function A(p, oj) from the imaginary time 
propagator using the methods of maximum entropy and singular value decomposition. We deter- 
mine the quasiparticle spectrum, which can be accurately parameterized by three parameters: an 
effective mass m* , a mean-field potential U, and a gap A. Large value of A/T c indicates that the 
system is not a BCS-type superfiuid at low temperatures. 

PACS numbers: 21.65.Cd, 03.75.Ss, 02.70.Ss, 26.60.Kp 



Dilute neutron matter is one of the simplest many- 
body nuclear systems. At sufficiently small densities its 
properties originate from the two-body s-wave interac- 
tion only. It is known that neutron matter has a posi- 
tive pressure at all densities (contrary to nuclear matter) 
which prevents fragmentation and it becomes superfiuid 
at low temperatures. From the theoretical point of view, 
pure and dilute neutron matter is a fascinating system 
since at a certain density range it becomes a nearly- 
universal Fermi gas. Such systems are presently of great 
interest as a result of an extraordinary progress in the 
field of cold atoms which have taken place over the last 
few years and in fact opened new chapter in many-body 
physics (see [l[ and references therein) . Taking advantage 
of the Feshbach resonances experimentalists can control 
the strength of the atom-atom interaction and achieve the 
so-called unitary regime. It corresponds to the situation 
where the average distance between fermionic atoms is 
large as compared to the interaction range ro, but much 
smaller than the scattering length a ie. pr 3 , <C 1 <C p\a\ 3 , 
where p is the particle number density. In the unitary 
regime the properties of dilute Fermi gases are universal, 
independent of the details of the interaction. Universal- 
ity of these systems make them fascinating theoretical 
playground, and obtained results turned out to be rel- 
evant to a wide range of fields like string theories, the 
quark-gluon plasma, and high T c superconductors. 

Since the Sq neutron- neutron interaction is charac- 
terized by the large scattering length a ss —18.5 fin, the 
unitary regime can be thought of as a limiting case of 
dilute neutron matter at the density range varying from 
0.001 to 0.01 fm~ 3 . One has to remember, however, that 
the influence of the effective range (r e ff « 2.8 fm)cannot 
be ignored since &F?~eff is of the order of unity |2j. The 
importance of other channels as well as of three-body 
forces is increasing with density. However at the density 
0.003 fm -3 , which we study in this paper, their influ- 
ence is marginal as compared to uncertainties of PIMC 



method and therefore will be neglected 0, 0] ■ 

Since even for the density p — 0.001 fm -3 dilute neu- 
tron matter is a strongly correlated Fermi gas (|fci?a| 3> 
1) only non-perturbative approaches are able to gain re- 
liable insight into physics of this system. The large class 
of such methods, which are known under the general 
name of Quantum Monte Carlo (QMC), have been used 
to date, although most of them concern the zero tem- 
perature properties [5j-[7|]. The finite temperature behav- 
ior has been studied in |8[ . This work presents the first 
ab initio, fully non-perturbative evaluation of thermal 
properties of low-density neutron matter (at about 2% of 
nuclear saturation density) free of uncontrolled approxi- 
mations within PIMC method. We focus on the effects 
generated by the finite effective range. 

Contrary to cold atomic gases, in order to capture 
physics of dilute neutron matter one has to use more re- 
alistic interaction than a simple contact, delta- like force. 
In the present paper we employ the two-body potential 
of the form: 



6g, r — r' = 
g, r - r' 6 Af b 
0. otherwise 



(1) 



where N b = {(±6,0,0), (0, ±6,0), (0,0, ±6)} represents 
the set of the nearest neighbor coordinates. This par- 
ticular form of the interaction is especially designed for 
the cubic lattice with the lattice constant b and enables 
to construct a fully non-perturbative approach without 
the sign problem (for more details see Ref. Q). It de- 
pends on two parameters (g and 6) which are adjusted to 
correctly reproduce the scattering length and the effec- 
tive range of neutron- neutron 1 So scattering amplitude 
[Io| . Hence we consider the system on a 3D spatial cubic 
lattice of length L = N s b with periodic boundary con- 
ditions. The lattice spacing b and size L introduce the 
natural ultraviolet (UV) and infrared (IR) momentum 
cut-offs given by p cu t = t/6 and po = 2ir / L, respectively. 
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The momentum space has the shape of a cubic lattice, 
with size 2ir/b and spacing 2n/L. To simplify the anal- 
ysis, however, we place the spherically symmetric UV 
cut-off, including momenta p < p cu t- 

To evaluate numerically expectation values of observ- 
ables we have followed the path integral approach de- 
scribed in Rcf. 
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Using Trotter expansion and subse- 
quently Hubbard-Stratonovich (H-S) transformation, the 
evaluation of the emerging path integral was performed 
using the Metropolis importance sampling. The crucial 
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consists m 



modification of the procedure described in 
the construction of such H-S transformation which allow 
to incorporate the off-site part of the interaction without 
generation of the sign problem. Namely, we have used 
the discrete H-S transformation of the form 



-tV 



n n 



i 



k 



(2) 



where <Ji are real numbers and h\(v) is the occupation 
number operator. The notable feature of this H-S trans- 
formation is the time reversal invariance of the corre- 
sponding imaginary time evolution operator. This prop- 
erty ensures that the probability measure used in the 
Metropolis algorithm is always positive 0, [l2[ . 

Calculations were performed on the lattice of size 
N s = 8 with the lattice constant b = 3.21 fm. The chem- 
ical potential was chosen in such a way to keep the to- 
tal number of particles between 53 and 57, which cor- 
responds to the density kp — 0.45 fm" 1 . The tempera- 
tures span the interval from 0.06 Ef (0.26 MeV) to I.Qef 
(4.3 MeV), where Ef is the Fermi energy. The number 
of imaginary time steps required to reach the conver- 
gence of the algorithm varies with temperature. At the 
lowest temperature 2360 imaginary time steps have been 
applied, whereas for the highest temperature only 216. 
The kinetic energy part of the Hamiltonian is defined in 
the restricted momentum space (p < p C ut) using the dis- 
persion relation of the form e(p) =p 2 /2m. Consequently 
during the imaginary time evolution the FFT algorithm 
has been used to switch between momentum and coordi- 



nate spaces [ll| . The number of generated uncorrected 
Monte Carlo samples allows to decrease the statistical 
error below 5%. At low temperatures the Singular Value 
Decomposition technique was applied to avoid instabili- 
ties of the algorithm. In all runs the single-particle oc- 
cupation probabilities for the highest energy states were 
below one percent at all temperatures. We have also 
performed a few exploratory simulations for the lattice 
of size N s = 10. The results were in a good agreement 
with those for N s = 8 lattice. 

In the Fig. [T]the low temperature behavior of the total 
energy and the chemical potential is presented for two 
different lattice sizes. The (shifted) total energy versus 
temperature for the free Fermi gas at the same particle 
density has also been plotted (solid line) . Note that after 
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FIG. 1: (Color online) Total energy E and chemical poten- 
tial )i as a function of temperature for dilute neutron matter 
at the density p = 0.003 fm" 3 (k F ~ 0.45 fm" 1 ). The total 
energy is denoted by red squares (8 3 lattice) and blue circles 
(10 3 lattice). The red up triangles and blue down triangles 
correspond to the chemical potential for 8 3 and 10 3 lattice, 
respectively. The solid line represents the energy of the nonin- 
teracting Fermi gas, shifted by a constant value. The dashed 
line shows an extrapolation of the energy and the chemical 
potential to T — limit. For comparison the total energy 
of the unitary Fermi gas is also plotted (open red squares). 
Dashed area for T = denotes the range where the results 
of other QMC results are located (see for example Ref. [(|). 
In the inset the rescaled condensate fraction is shown as a 
function of temperature, red squares and blue circles denote 
8 3 and 10 3 lattices, respectively. Crosspoint determines the 
critical temperature of the superfluid-normal phase transition, 
T c w 0.09 ef- 



shifting of the free Fermi gas energy by 0.52_Effg the 
curve reproduces Monte Carlo results for T > 0.15ef 
(£"ffg = §A^£j? is the free Fermi gas energy at T = 0). 
Below this temperature the deviation from the free Fermi 
gas behavior is clearly visible. The chemical potential is 
approximately constant for T < 0.1 ef- 

The critical temperature of the superfluid-normal 
phase transition has been determined using the method 
based on the finite size scaling of the correlation func- 
tion. Similar technique was used to determine the critical 
temperature at the unitary limit (see Refs. 



11, 131 for 



details). The volume-dependent estimation of the critical 
temperature Tc was obtained by finding the crossing 
point of the rescaled condensate fraction for two different 
lattice sizes Nij. As N^j — > oo, the series converges 
to T c and one can extract the limiting value. We have 
determined T c using results for two lattices Nij = 8, 10. 
Such large lattices and rather small filling factor which 
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in both cases reads v — N/2N^ « 5% are enough to es- 
timate the critical temperature with uncertainty smaller 
than 20% (in fact this procedure applied to the unitary 
gas gives estimation of the T c with the relative error 
smaller than 10%). The estimate of the critical temper- 
ature reads T c k, 0.09 ep. Note that T c is considerably 
lower than the temperature for the onset of deviation 
from the free Fermi gas behavior. 

Within the PIMC framework one cannot reach directly 
the T = limit. However the ground state energy can 
be obtained by performing an extrapolation of results 
to zero temperature limit. In our case this procedure 
provides the ground state energy E/E-pfg — 0.46(2) 
(E/N = 1.22(5) MeV). This value is considerably lower 
(by about 20%-40%) than values obtained by other MC 
calculations (see for example Ref. @). This is most likely 
due to the fact, that our approach is based on fully unre- 
stricted path integral calculations and, within statistical 
errors due to the Monte Carlo procedure, gives essentially 
exact results. 

The gap in the fermionic spectrum, related to superflu- 
idity, has been computed from the spectral weight func- 
tion A(p, ui) by performing the analytic continuation of 
the imaginary time propagator Q(p, r) to real frequencies 
jl4] . This procedure is equivalent to solving the integral 
equation: 



a(p,r) = -i- 

Z7T 



duiA(p, ui) 



exp(— cut) 
1 + exp(— w/3) ; 



(3) 



where Q(p,r) is known from the Monte Carlo calcula- 
tions for 51 different values of r € [Q,/3 = 1/T]. The in- 
verse problem is however numerically ill-posed i.e. there 
is an infinite class of solutions for A(p,u>) which satisfy 
Eq. ([3]) within uncertainties generated by the Monte 
Carlo method. Therefore we have used two independent 
methods based on completely different mathematical ap- 
proaches. 

The first one, the maximum entropy method, is based 
on Bayes' theorem [lg. It treats the values of G(p,Ti) 
(i = 0,1, ...,50) provided by QMC simulation as nor- 
mally distributed random numbers, around the true val- 
ues Q(p,Ti), and searches for the most probable solution 
assuming some a priori knowledge concerning the spec- 
tral function. As an a priori information we have used 
constraints: 

/+oo I 
^A(p,u,) = l, (4) 
-oo 2lT 

/ — A(p,w)— r — = n{p), 5 

and we have assumed a Gaussian-like structure for 
A(p,u). In the formula ([5]) n(p) represents the occu- 
pation probability of the state with momentum p which 
is known from the Monte Carlo simulation. 
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FIG. 2: (Color online) Spectral weight function A(p, ui) at the 
temperature T = 0.06 ef and lattice size N 3 — 10 obtained 
by the maximum entropy method. Points indicate localiza- 
tions of maxima for fixed values of momenta. Dashed lines 
correspond to the fit of the BCS-type formula given by ([6]). 
In the inset the spectral weight function at the Fermi level is 
presented. 



The second method is based on the singular value de- 
composition (SVD) of the integral kernel K, of Eq. ([3]), 
which can be rewritten in the operator form as G(p, Tj) = 
(fCA)(p,Ti), The operator K, possesses the singular sys- 
tem which forms a suitable basis for the expansion of the 
projected spectral weight function A(p,u>) onto a sub- 
space where the inverse problem is well-posed (l6j . Since 
the method provides only projection of the "true" solu- 
tion, it does not require any a priori information, con- 
trary to the maximum entropy method. However, since 
Q(p, n) include statistical errors due to the Monte Carlo 
procedure, the projected solution A(p,uj) is also affected 
by this uncertainty. One can use this flexibility by choos- 
ing the solution satisfying the constraints (j4|) [rjj. The 
details of both methods will be discussed elsewhere [3| • 

The spectral weight function for the lowest tempera- 
ture T — 0.06 Ef obtained for A^ s = 10 lattice is shown 
in the Fig. [2] The same outcome has been generated 
by both methods (maximum entropy and SVD) indepen- 
dently. The presence of a "pairing" gap is clearly visible 
for this temperature. 

Figure [3] presents the quasiparticle excitation spec- 
trum extracted from the spectral weight function for 
T = 0.06 Ef- We have found that the quasiparticle exci- 
tations can be accurately parameterized by the BCS-like 
formula: 



E(p) = ±y 



2m* 



U 



(6) 



where m* is an effective mass, U the mean field potential 
and A is the "pairing" gap. The values of these parame- 
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FIG. 3: (Color online) Quasiparticle energies E(p) (squares) 
extracted from the spectral weight function A(p, w) at T = 
0.06 ef- The line denotes results obtained under assumption 
that the system is composed of independent quasiparticles. 
Dashed line corresponds to quasiparticle energies at the uni- 
tary limit. 



ters were estimated as m* jm = 1.1(1), U/ef = —0.26(6) 
and A/e F = 0.25(5). 

Note that the ratio A/T c « 2.8 is significantly higher 
than the well-known value 1.76 predicted by BCS theory. 
The similar deviation from the BCS value is typical for 
high-temperature superconductors [19|_and also for cold 
atomic gases in the unitary regime [ll[. Therefore we 
conclude, that the dilute neutron matter at this density 
is not a BCS-type superfluid. Note also that to estimate 
the value of A/T c we have used the value of the energy 
gap at the temperature T = 0.06 ep, which is expected 
to be slightly lower than the value of the gap at zero 
temperature. 

It is instructive to compare quasiparticle excitation en- 
ergies with those extracted from the susceptibility func- 
tion under assumption that the system is composed of 
independent quasiparticles. Under this assumption the 
imaginary time propagator is simply given by: 



G(p,r) = - 



= -tS(p) 



1 



and one can easily evaluate the susceptibility: 

rP i e /3£(p) . 

X(P) = / drg{p,T) = 



E(p) e^P) 



1 



(7) 



(8) 



From the calculated one-body propagator within the 
Monte Carlo algorithm one can extract the spectrum of 
the elementary fermionic excitations inverting the Eq. 
(J5J). The extracted spectrum of quasiparticle energies 
turns out to reproduce very well (within error bars) the 
quasiparticle spectrum derived from the spectral func- 
tion, see Fig. [3l The same property is shared by unitary 
cold atomic gas at temperatures below the critical tem- 
perature [20l | . 



Comparison of our results with those obtained in the 
limit r e ff — > provides an information about the influence 



of the effective range. From the data reported in Ref. [11 1 
we infer that the effects of the effective range do not sig- 
nificantly alter the ground state energy. The value of the 
energy gap and the critical temperature decreases consid- 
erably (at r cff -> 0: A^/e F « 0.41 and T c (0) /£f ps 0.13). 
However, surprising ly the ratio /T c (0) « 3.2 remains 
approximately constant (taking into account uncertain- 
ties of our estimation) when increasing r e ff to the value 
associated with 1 Sq neutron-neutron interaction. Note 
also that the equation of state exhibits the existence of 
the second temperature scale, which can be attributed 
to the onset of deviations of E/Eyfg from the (shifted) 
energy of the free Fermi gas. It bears similarity to the 
case of the unitary Fermi gas, where the existence of the 
so-called "pseudogap" above T c is reported [2fj| . 

Summarizing, our results do not indicate the presence 
of qualitative changes in comparison to the case of zero 
effective range. In conclusion the main aspects of physics 
at the unitary regime survive in the limit of dilute neu- 
tron matter. 
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